FunTemp_d0  =   @(x)x.^0.5;
FunTemp_d1  =   @(x)0.5*x.^(-0.5);
FunTemp_d2  =   @(x)-0.25*x.^(-1.5);

MidLeft     =   20;
MidRight    =   90;
x_left      =   linspace(0,MidLeft,3)';
x_middle    =   linspace(MidLeft+1,MidRight-1,30)';
x_right     =   linspace(MidRight,100,3)';
x           =   [x_left;x_middle;x_right];

s_left      =   2;
yd0_left    =   x_left*s_left;
yd1_left    =   ones(size(yd0_left))*s_left;
yd2_left    =   zeros(size(yd0_left));

LinCoef     =   s_left-FunTemp_d1(MidLeft);
Const       =   yd0_left(end)-FunTemp_d0(MidLeft);
yd0_middle  =   Const+...
                LinCoef*(x_middle-MidLeft)+FunTemp_d0(x_middle);
yd1_middle  =   LinCoef+FunTemp_d1(x_middle);
yd2_middle  =   FunTemp_d2(x_middle);

s_right     =   LinCoef+FunTemp_d1(MidRight);   
yd0_right   =   Const+...
                LinCoef*(MidRight-MidLeft)+FunTemp_d0(MidRight)+...
                s_right*(x_right-MidRight);       
yd1_right   =   ones(size(yd0_right))*s_right;
yd2_right   =   zeros(size(yd0_right));

yd0         =   [yd0_left;yd0_middle;yd0_right];
yd1         =   [yd1_left;yd1_middle;yd1_right];
yd2         =   [yd2_left;yd2_middle;yd2_right];

Breaks      =   [linspace(1,MidLeft,2)';...
                 linspace(MidLeft,MidRight,10)';...
                 linspace(MidRight,99,2)'];
k           =   3;
Node        =   Spline1dNode(Breaks,k);
Basis       =   Spline1dBasis(Breaks,k,Node,[0;1;2]);

Yd0_Node    =   Node;
Yd0_Node(Node<=MidLeft)     =   s_left*Node(Node<=MidLeft);
Yd0_Node(Node>MidLeft & Node<MidRight)  =   ...
                Const+...
                LinCoef*(Node(Node>MidLeft & Node<MidRight)-MidLeft)+...
                FunTemp_d0(Node(Node>MidLeft & Node<MidRight));
Yd0_Node(Node>=MidRight)    =   ...
                Const+...
                LinCoef*(MidRight-MidLeft)+FunTemp_d0(MidRight)+...
                s_right*(Node(Node>=MidRight)-MidRight);  
            
BCoef       =   Basis{1}\Yd0_Node;         
[PC0,~,UniBreaks]       =   BCoef2PCoef(Breaks,k,BCoef);
PC1         =   PCoefDer(PC0,1);
PC2         =   PCoefDer(PC1,1);

xx          =   linspace(x(1),x(end),100)';
Bxx         =   Spline1dBasis(Breaks,k,xx,[0;1;2]);

figure;
subplot(1,3,1)
plot(x,yd0,'k-');
hold on
plot(xx,Bxx{1}*BCoef,'ro');
subplot(1,3,2)
plot(x,yd1,'k-');
hold on
plot(xx,Bxx{2}*BCoef,'ro');
subplot(1,3,3)
plot(x,yd2,'k-');
hold on
plot(xx,Bxx{3}*BCoef,'ro');
figure;
subplot(1,3,1)
plot(x,yd0,'k-');
hold on
plot(xx,Spline1dPPEval(UniBreaks,PC0,xx,0),'ro');
subplot(1,3,2)
plot(x,yd1,'k-');
hold on
plot(xx,Spline1dPPEval(UniBreaks,PC1,xx,0),'ro');
subplot(1,3,3)
plot(x,yd2,'k-');
hold on
plot(xx,Spline1dPPEval(UniBreaks,PC2,xx,0),'ro');